function [phi,psi,x,J,er] = sn_one_d(in)
% function [phi,psi,x,J] = sn_one_d(in)
%
% A simple 1-D discrete ordinates solver to test acceleration schemes
% including 
%   1. coarse mesh rebalance
%   2. coarse mesh finite difference
%   3. a new(?) approach based on sparse characteristics in homogenized
%      coarse regions

numg = in.numg;
ord  = in.ord;
Q = zeros(  sum(in.xfm), ord, numg );

bound       =   [ord/2+1 ord; 1 ord/2];
gbound      =   [ 1  numg];
git         =   1;

% -------------------------------------------------------------------------
% -------------------------------------- Discretizations
[mu,w] = S_1D(in.ord);
j = 0;
for i = 1:length(in.xfm)
    dx( (j+1):(j+in.xfm(i))   )  = (in.xcm(i+1)-in.xcm(i))/in.xfm(i);
    for g=gbound(1):git:gbound(2)
        % src is the volumetric angular source, is n/cm^3-s
        % so we need to divide by 2 to distribute in angle
        Q( (j+1):(j+in.xfm(i)), :, g)  = in.src(g,i)/2;
    end
    mtt( (j+1):(j+in.xfm(i))   )  = in.mt(i);  % assign mat to each f mesh
    j = sum(in.xfm(1:i));
end
n = sum(in.xfm);
xx = zeros(length(dx)+1);
for i = 2:length(xx)
    xx(i) = xx(i-1) + dx(i-1);  % actual x's
end

x(1)=dx(1)/2;
for i = 2:length(dx)
    x(i) = x(i-1) + 0.5*(dx(i-1)+dx(i));
end
xc = 0.5*(in.xcm(2:end)+in.xcm(1:end-1));

% -------------------------------------------------------------------------
% -------------------------------------- Matrix Pre-allocation
psi         =   zeros(n+1,ord,numg);        % Angular Flux
%cur         =   zeros(n,numg);              % Current
phi_o       =   zeros(n,numg);              % Scalar Flux Old
phi         =   zeros(n,numg);              % Scalar Flux New
con1        =   zeros(n,ord,numg);          % Transport Coef.
con2        =   zeros(n,ord,numg);          % Transport Coef.
con3        =   zeros(n,ord,numg);          % Transport Coef.
%q           =   zeros(n,ord,numg);          % Scattered Source

% ------------- Precompute Coefficients
for g = gbound(1):git:gbound(2)
    for k = 1:n
        m = mtt(k);
        con1(k,:,g) = in.data((m-1)*numg+g,1)*dx(k)./(2.0*mu(:));
        con2(k,:,g) = (1.0-sign(mu(:)').*con1(k,:,g))./...
            (1.0+sign(mu(:)').*con1(k,:,g));
        con3(k,:,g) = sign(mu(:)').*dx(k)./(mu(:)'.*...
            (1.0+sign(mu(:)').*con1(k,:,g)));
    end
end

q = Q;

% ----------------- Solution Algorithm ------------------------------------
tic
for g = gbound(1):git:gbound(2)
    
    % ----------------- Convergence Parameters ----------------------------
    err  = 10;
    iter = 0;
    
    while err > in.maxerr && iter < in.maxit
        
        % reflective condition
        if in.bcL == 1
            psi(1,bound(1,1):bound(1,2),:)=...
                psi(1,bound(2,2):-1:bound(2,1),:);
        end
        if in.bcR == 1
            psi(n+1,bound(2,2):-1:bound(2,1),:)=...
                psi(n+1,bound(1,1):bound(1,2),:);
        end
        % -------------------------------------- Mu > 0
        for k = 1:1:n
            for j = bound(1,1):bound(1,2)
                psi(k+1,j,g) = psi(k,j,g)*con2(k,j,g)+q(k,j,g)*con3(k,j,g);
            end
        end
        % -------------------------------------- Mu < 0
        for k = n:-1:1
            for j = bound(2,1):bound(2,2)
                psi(k,j,g) = psi(k+1,j,g)*con2(k,j,g)+q(k,j,g)*con3(k,j,g);
            end
        end
        
        % ---- Scalar Flux (use average of edge angular fluxes)
        for k = 1:n
            phi(k,g) = sum(0.5*w'.*(psi(k,:,g)+psi(k+1,:,g)));
        end
        
%         plot(x,phi,'k')
%         grid on
        % acceleration schemes (update both psi and phi)
        if ( in.acc == 1 )
            [tmp,phi] = cmr(psi,phi,Q,in,mu,w,dx);
        elseif ( in.acc == 2 )
            [tmp,phi] = cmfd(psi,phi,Q,in,mu,w,dx);
        elseif ( in.acc == 3)
            [tmp,phi] = cmsc(psi,phi,Q,in,mu,w,dx);
        end
%         hold on
%         plot(x,phi,'b--')
%         lala=1;
        % ---- Updated Source Term
        for z   =   g:git:gbound(2) % only down scattering
            for kk = 1:n
                q(kk,:,z) = Q(kk,:,z); % reset
            end
            if in.adj == 0
                for k = 1:n
                    m = mtt(k);
                    for gg = gbound(1):git:gbound(2) % group gg to group z
                        q(k,:,z) = q(k,:,z) + ...
                            0.5*in.data((m-1)*numg+gg,2+z)*phi(k,gg);
                    end
                end
            else
                for k = 1:n
                    m = mtt(k);
                    for gg = gbound(1):git:gbound(2) % group gg to group z
                        q(k,:,z) = q(k,:,z) + ...
                            0.5*in.data((m-1)*numg+z,2+gg)*phi(k,gg);
                    end
                end
            end
        end
        % ---- Scalar Flux Error Between Iterations
        err = max(max(abs(phi-phi_o)./phi));
        % ---- Reset Scalar Flux
        phi_o = phi;
        % ---- Iteration Counter
        iter = iter + 1;
        if (mod(iter,100)==0)
            disp(['   iter = ',num2str(iter),'   err_phi = ',num2str(err)])
        end
       % disp(['   iter = ',num2str(iter),'   err_phi = ',num2str(err)])
        er(iter)=err;

    end
    
end

% current density
J = zeros(n,in.numg);
for k = 1:n
    J(k,g) = sum(w'.*mu'*0.25.*(psi(k,:,g)+psi(k+1,:,g)));
end

% % Some stuff for plotting
% x(1)=dx(1)/2;
% for i = 2:length(dx)
%     x(i) = x(i-1) + 0.5*(dx(i-1)+dx(i));
% end
% if in.acc == 0
%     figure(1)
%     plot(x,phi,'k',x,J,'k','LineWidth',2)
% elseif in.acc ==1
%     figure(1)
%     plot(x,phi,'g--',x,J,'g--','LineWidth',2)
% else
%     figure(1)
%     plot(x,phi,'r.',x,J,'r.','LineWidth',2)    
% end

disp(['pi iterations: ',num2str(iter)])

end

